Cargo librerias

library(spatstat.core)
## Loading required package: spatstat.data
## Loading required package: spatstat.geom
## spatstat.geom 2.3-0
## Loading required package: nlme
## Loading required package: rpart
## spatstat.core 2.3-1
library(plotly)
## Loading required package: ggplot2
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout

Principio de Aceptacion-Rechazo para un Proceso Binomial

# Graficacion
cuadri <- function(x, y){sqrt(x^2+y^2)/sqrt(8)} # Funcion de intensidad, toma valores entre 0 y 1
n<-10 # cant de datos
x <- y <- seq(0, 2, length= n)
xx<-rep(x,n)
yy<-rep(y,rep(n,n))
cbind(xx,yy) # grilla en el soporte
##               xx        yy
##   [1,] 0.0000000 0.0000000
##   [2,] 0.2222222 0.0000000
##   [3,] 0.4444444 0.0000000
##   [4,] 0.6666667 0.0000000
##   [5,] 0.8888889 0.0000000
##   [6,] 1.1111111 0.0000000
##   [7,] 1.3333333 0.0000000
##   [8,] 1.5555556 0.0000000
##   [9,] 1.7777778 0.0000000
##  [10,] 2.0000000 0.0000000
##  [11,] 0.0000000 0.2222222
##  [12,] 0.2222222 0.2222222
##  [13,] 0.4444444 0.2222222
##  [14,] 0.6666667 0.2222222
##  [15,] 0.8888889 0.2222222
##  [16,] 1.1111111 0.2222222
##  [17,] 1.3333333 0.2222222
##  [18,] 1.5555556 0.2222222
##  [19,] 1.7777778 0.2222222
##  [20,] 2.0000000 0.2222222
##  [21,] 0.0000000 0.4444444
##  [22,] 0.2222222 0.4444444
##  [23,] 0.4444444 0.4444444
##  [24,] 0.6666667 0.4444444
##  [25,] 0.8888889 0.4444444
##  [26,] 1.1111111 0.4444444
##  [27,] 1.3333333 0.4444444
##  [28,] 1.5555556 0.4444444
##  [29,] 1.7777778 0.4444444
##  [30,] 2.0000000 0.4444444
##  [31,] 0.0000000 0.6666667
##  [32,] 0.2222222 0.6666667
##  [33,] 0.4444444 0.6666667
##  [34,] 0.6666667 0.6666667
##  [35,] 0.8888889 0.6666667
##  [36,] 1.1111111 0.6666667
##  [37,] 1.3333333 0.6666667
##  [38,] 1.5555556 0.6666667
##  [39,] 1.7777778 0.6666667
##  [40,] 2.0000000 0.6666667
##  [41,] 0.0000000 0.8888889
##  [42,] 0.2222222 0.8888889
##  [43,] 0.4444444 0.8888889
##  [44,] 0.6666667 0.8888889
##  [45,] 0.8888889 0.8888889
##  [46,] 1.1111111 0.8888889
##  [47,] 1.3333333 0.8888889
##  [48,] 1.5555556 0.8888889
##  [49,] 1.7777778 0.8888889
##  [50,] 2.0000000 0.8888889
##  [51,] 0.0000000 1.1111111
##  [52,] 0.2222222 1.1111111
##  [53,] 0.4444444 1.1111111
##  [54,] 0.6666667 1.1111111
##  [55,] 0.8888889 1.1111111
##  [56,] 1.1111111 1.1111111
##  [57,] 1.3333333 1.1111111
##  [58,] 1.5555556 1.1111111
##  [59,] 1.7777778 1.1111111
##  [60,] 2.0000000 1.1111111
##  [61,] 0.0000000 1.3333333
##  [62,] 0.2222222 1.3333333
##  [63,] 0.4444444 1.3333333
##  [64,] 0.6666667 1.3333333
##  [65,] 0.8888889 1.3333333
##  [66,] 1.1111111 1.3333333
##  [67,] 1.3333333 1.3333333
##  [68,] 1.5555556 1.3333333
##  [69,] 1.7777778 1.3333333
##  [70,] 2.0000000 1.3333333
##  [71,] 0.0000000 1.5555556
##  [72,] 0.2222222 1.5555556
##  [73,] 0.4444444 1.5555556
##  [74,] 0.6666667 1.5555556
##  [75,] 0.8888889 1.5555556
##  [76,] 1.1111111 1.5555556
##  [77,] 1.3333333 1.5555556
##  [78,] 1.5555556 1.5555556
##  [79,] 1.7777778 1.5555556
##  [80,] 2.0000000 1.5555556
##  [81,] 0.0000000 1.7777778
##  [82,] 0.2222222 1.7777778
##  [83,] 0.4444444 1.7777778
##  [84,] 0.6666667 1.7777778
##  [85,] 0.8888889 1.7777778
##  [86,] 1.1111111 1.7777778
##  [87,] 1.3333333 1.7777778
##  [88,] 1.5555556 1.7777778
##  [89,] 1.7777778 1.7777778
##  [90,] 2.0000000 1.7777778
##  [91,] 0.0000000 2.0000000
##  [92,] 0.2222222 2.0000000
##  [93,] 0.4444444 2.0000000
##  [94,] 0.6666667 2.0000000
##  [95,] 0.8888889 2.0000000
##  [96,] 1.1111111 2.0000000
##  [97,] 1.3333333 2.0000000
##  [98,] 1.5555556 2.0000000
##  [99,] 1.7777778 2.0000000
## [100,] 2.0000000 2.0000000
z <- outer(x, y, cuadri) # intensidad calculada en cada punto de grilla
u<-runif(n*n) # puntos aleatorios x y
graf<-persp(x, y, z,theta = 65,phi=25,zlab="u")
cumplen<-(u<z) # restricción para simular bajo la intensidad
points(trans3d(xx, yy, u,pmat=graf), col = 2+cumplen, pch = 16)

#  Simulacion
set.seed(1) # fijo la semilla
cant<-1000 # cant de obs
x<-runif(cant,0,2) # genero x entre 0 y 2
y<-runif(cant,0,2) # genero y entre 0 y 2
u<-runif(cant,0,1) # generu u entre 0 y 1
z<-cuadri(x,y) # evaluo en la funcion de intensidad
cumplen<-(u<z) # Restriccion
sum(cumplen) # cuantos quedan
## [1] 554
# Grafico 100 casos
par(pty="s")
plot(x[cumplen][1:100],y[cumplen][1:100],xlab = "",ylab="",xaxt='n',yaxt='n')

# Grafico todos los casos
par(pty="s")
plot(x[cumplen],y[cumplen],xlab = "",ylab="",xaxt='n',yaxt='n')

Simulacion de Complete Spatial Randomness (CSR) o Independent Random Process (IRP) Funciones de spatstat,core (Uniform) Binomial Point Process

set.seed(1)
cant<-100
sim.BPP<-runifpoint(cant)
class(sim.BPP)
## [1] "ppp"
plot(sim.BPP)

set.seed(2)
sim.BPP<-runifpoint(cant)
plot(sim.BPP)

# muchas realizaciones
sim.BPP<-runifpoint(cant,nsim = 6)
plot(sim.BPP)

class(sim.BPP)
## [1] "ppplist" "solist"  "anylist" "listof"  "list"

Proceso de poisson homogeneo

set.seed(1)
lambda<-100
sim.CSR<-rpoispp(lambda)
plot(sim.CSR,main = paste("Poisson Homogeneo lambda =",lambda))

set.seed(2)
lambda<-100
sim.CSR<-rpoispp(lambda)
plot(sim.CSR,main = paste("Poisson Homogeneo lambda =",lambda))

Proceso de poisson INhomogeneo

set.seed(1)
sim.INHa <- rpoispp(function(x,y) {800*(x-0.5)^2 + 800*(y-0.5)^2},nsim=1)
plot(sim.INHa)

set.seed(1)
sim.INHb <- rpoispp(function(x,y) {200-800*(x-0.5)^2 + 200-800*(y-0.5)^2},nsim=1)
plot(sim.INHb)

Procesos de Superposición (Poisson)

# a lo bruto
plot(sim.INHa,type="n",main="Superposición")
points(sim.INHa,col="blue")
points(sim.INHb,col="red")

# mas sofisticado
super<-superimpose(sim.INHa,sim.INHb) # superposicion de ppps
plot(super)

Analisis de Concentracion PP Homogeneo

# El proceso CSR
sim.CSR
## Planar point pattern: 91 points
## window: rectangle = [0, 1] x [0, 1] units
plot(sim.CSR)

# Quadrat Count
Q <- quadratcount(sim.CSR, nx= 4, ny=4)
plot(sim.CSR, pch=20, cols="grey70", main=NULL)  # Plot points
plot(Q, add=TRUE)  # Add quadrat grid

# Quadrat Density
Q.d <- intensity(Q)
plot(intensity(Q, image=TRUE), main=NULL, las=1)  # Plot density raster
plot(sim.CSR, pch=1, cex=0.01, col="grey", add=TRUE)  # Add points

Test de uniformidad de distribucion PP Homogeneo

QT <- quadrat.test(sim.CSR, nx=4, ny=4)
QT
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  sim.CSR
## X2 = 16.077, df = 15, p-value = 0.7539
## alternative hypothesis: two.sided
## 
## Quadrats: 4 by 4 grid of tiles
plot(QT)

Analisis de Concentracion PP InHomogeneo

# Quadrat Count
Q <- quadratcount(sim.INHa, nx= 4, ny=4)
plot(sim.INHa, pch=20, cols="grey70", main=NULL)  # Plot points
plot(Q, add=TRUE)  # Add quadrat grid

# Quadrat Density
Q.d <- intensity(Q)
plot(intensity(Q, image=TRUE), main=NULL, las=1)  # Plot density raster
plot(sim.INHa, pch=1, cex=0.01, col="grey", add=TRUE)  # Add points

Test de uniformidad de distribucion PP InHomogeneo

QT <- quadrat.test(sim.INHa, nx=4, ny=4)
QT
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  sim.INHa
## X2 = 73.619, df = 15, p-value = 2.009e-09
## alternative hypothesis: two.sided
## 
## Quadrats: 4 by 4 grid of tiles
plot(QT)

Proceso de Thomas

par(mfrow=c(2,2))
par(mar=c(1,1,1,1))
set.seed(1)
inten<-30
disper<-0.02
mu<-10
sim.THO <- rThomas(inten,disper,mu)
plot(sim.THO,main=paste("Lambda=",inten,"Mu=",mu,"Sigma=",disper))
set.seed(1)
inten<-30
disper<-0.05
mu<-10
sim.THO <- rThomas(inten,disper,mu)
plot(sim.THO,main=paste("Lambda=",inten,"Mu=",mu,"Sigma=",disper))
set.seed(1)
inten<-10
disper<-0.02
mu<-30
sim.THO <- rThomas(inten,disper,mu)
plot(sim.THO,main=paste("Lambda=",inten,"Mu=",mu,"Sigma=",disper))
set.seed(1)
inten<-10
disper<-0.05
mu<-30
sim.THO <- rThomas(inten,disper,mu)
plot(sim.THO,main=paste("Lambda=",inten,"Mu=",mu,"Sigma=",disper))

par(mfrow=c(1,1))

Proceso de Secuencial Inhibicion Simple (SSI)

par(mfrow=c(2,2))
par(mar=c(1,1,1,1))
set.seed(1)
delta<-0.01
sim.SSI <- rSSI(delta,100)
plot(sim.SSI,main=paste("n=100 delta=",delta))
set.seed(1)
delta<-0.05
sim.SSIa <- rSSI(delta,100)
plot(sim.SSIa,main=paste("n=100 delta=",delta))
set.seed(1)
delta<-0.1
sim.SSI <- rSSI(delta,100)
## Warning in rSSI(delta, 100): Gave up after 1000 attempts with only 69 points
## placed out of 100
plot(sim.SSI,main=paste("n=100 delta=",delta))
set.seed(1)
delta<-0.2
sim.SSIb <- rSSI(delta,100)
## Warning in rSSI(delta, 100): Window is too small to fit 100 points at minimum
## separation 0.2 (absolute maximum number is 31)
## Warning in rSSI(delta, 100): Gave up after 1000 attempts with only 21 points
## placed out of 100
plot(sim.SSIb,main=paste("n=100 delta=",delta))

Procesos de COX

library(RandomFields)
## Loading required package: sp
## Loading required package: RandomFieldsUtils
## 
## Attaching package: 'RandomFields'
## The following object is masked from 'package:RandomFieldsUtils':
## 
##     RFoptions
## The following object is masked from 'package:nlme':
## 
##     Variogram
library(RandomFieldsUtils)
# Defino el campo medio de la funcion de intensidad
m <- as.im(function(x, y){5 - 1.5 * (x - 0.5)^2 + 2 * (y - 0.5)^2}, W=owin())
class(m)
## [1] "im"
plot(m)

# Genero un Proceso de Cox Log Gausiano a partir del campo medio generado previamente
X <- rLGCP("gauss", m, var=0.15, scale =0.5)
# Grafico del campo aleatorio + la realizacion del proceso
plot(attr(X, "Lambda"))
points(X)

# Probemos tra ventana
#
ho <- owin(poly=list(list(x=c(0,2,1,0), y=c(0,0.3,1,0.8)),
                       list(x=c(0.6,0.4,0.4,0.6), y=c(0.2,0.2,0.4,0.4))))

m <- as.im(function(x, y){5 - 1.5 * (x - 0.5)^2 + 2 * (y - 0.5)^2}, W=ho)
plot(m)

La Funcion G y F en CSR

par(mfrow=c(1,3))
plot(sim.CSR)
env.CSR <- envelope(sim.CSR, fun = Gest, nsim = 200,nrank=5)
## Generating 200 simulations of CSR  ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198. 200.
## 
## Done.
plot(env.CSR)
env.CSR <- envelope(sim.CSR, fun = Fest, nsim = 200,nrank=5)
## Generating 200 simulations of CSR  ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198. 200.
## 
## Done.
plot(env.CSR)

par(mfrow=c(1,1))

La Funcion G y F en Thomas (Concentracion)

par(mfrow=c(1,3))
plot(sim.THO)
env.THO <- envelope(sim.THO, fun = Gest, nsim = 200,nrank=5)
## Generating 200 simulations of CSR  ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198. 200.
## 
## Done.
plot(env.THO)
env.THO <- envelope(sim.THO, fun = Fest, nsim = 200,nrank=5)
## Generating 200 simulations of CSR  ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198. 200.
## 
## Done.
plot(env.THO)

par(mfrow=c(1,1))

La Funcion G y F en SSI (Regularidad)

par(mfrow=c(1,3))
sim.SSI<-sim.SSIa
plot(sim.SSI)
env.SSI <- envelope(sim.SSI, fun = Gest, nsim = 200,nrank=5)
## Generating 200 simulations of CSR  ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198. 200.
## 
## Done.
plot(env.SSI)
env.SSI <- envelope(sim.SSI, fun = Fest, nsim = 200,nrank=5)
## Generating 200 simulations of CSR  ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198. 200.
## 
## Done.
plot(env.SSI)

par(mfrow=c(1,1))

La Funcion G y F en CSR

par(mfrow=c(1,3))
plot(sim.CSR)
env.CSR <- envelope(sim.CSR, fun = Kest, nsim = 200,nrank=5)
## Generating 200 simulations of CSR  ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198. 200.
## 
## Done.
plot(env.CSR)
env.CSR <- envelope(sim.CSR, fun = Lest, nsim = 200,nrank=5)
## Generating 200 simulations of CSR  ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198. 200.
## 
## Done.
plot(env.CSR, . -r ~ r)

par(mfrow=c(1,1))

La Funcion K y L en Thomas

par(mfrow=c(1,3))
plot(sim.THO)
env.THO <- envelope(sim.THO, fun = Kest, nsim = 200,nrank=5)
## Generating 200 simulations of CSR  ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198. 200.
## 
## Done.
plot(env.THO)
env.THO <- envelope(sim.THO, fun = Lest, nsim = 200,nrank=5)
## Generating 200 simulations of CSR  ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198. 200.
## 
## Done.
plot(env.THO, . -r ~ r)

par(mfrow=c(1,1))

La Funcion K y L en SSI

par(mfrow=c(1,3))
sim.SSI<-sim.SSIa
plot(sim.SSI)
env.SSI <- envelope(sim.SSI, fun = Kest, nsim = 200,nrank=5)
## Generating 200 simulations of CSR  ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198. 200.
## 
## Done.
plot(env.SSI)
env.SSI <- envelope(sim.SSI, fun = Lest, nsim = 200,nrank=5)
## Generating 200 simulations of CSR  ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198. 200.
## 
## Done.
plot(env.SSI, . -r ~ r)

par(mfrow=c(1,1))

Fry Plots o graficos de Petterson

par(mfrow=c(2,3))
par(mar=c(1,1,1,1))
plot(sim.CSR)
plot(sim.THO)
plot(sim.SSIb)
fryplot(sim.CSR)
fryplot(sim.THO)
fryplot(sim.SSIb)

par(mfrow=c(1,1))

Estimcion de la intensidad para poisson con modelo parametrico (funcion ppm)

Caso Homogeneo

# modelo homogeneo
ajus.CSR.1<-ppm(sim.CSR~1)
ajus.CSR.1
## Stationary Poisson process
## Intensity: 91
##             Estimate      S.E.  CI95.lo CI95.hi Ztest     Zval
## log(lambda)  4.51086 0.1048285 4.305399 4.71632   *** 43.03086
# modelo cinhomogeneo on tendencia espacial lineal
ajus.CSR.2<-ppm(sim.CSR~x+y)
ajus.CSR.2
## Nonstationary Poisson process
## 
## Log intensity:  ~x + y
## 
## Fitted trend coefficients:
## (Intercept)           x           y 
##  4.52381507  0.06116668 -0.08803713 
## 
##                Estimate      S.E.    CI95.lo   CI95.hi Ztest       Zval
## (Intercept)  4.52381507 0.2769819  3.9809404 5.0666897   *** 16.3325273
## x            0.06116668 0.3633154 -0.6509185 0.7732518        0.1683569
## y           -0.08803713 0.3633713 -0.8002318 0.6241576       -0.2422787
plot(ajus.CSR.2,pause=FALSE,superimpose = FALSE)

plot(ajus.CSR.2,pause=FALSE,how="persp", theta=-30,phi=40,d=4)

plot(ajus.CSR.2,pause=FALSE,how="contour", theta=-30,phi=40,d=4)

pred2<-predict(ajus.CSR.2)
plot(pred2)

Caso InHomogeneo

# modelo homogeneo
ajus.INH.1<-ppm(sim.INHa~1)
ajus.INH.1
## Stationary Poisson process
## Intensity: 126
##             Estimate       S.E.  CI95.lo  CI95.hi Ztest     Zval
## log(lambda) 4.836282 0.08908708 4.661674 5.010889   *** 54.28713
# modelo cinhomogeneo on tendencia espacial lineal
ajus.INH.2<-ppm(sim.INHa~x+y)
ajus.INH.2
## Nonstationary Poisson process
## 
## Log intensity:  ~x + y
## 
## Fitted trend coefficients:
## (Intercept)           x           y 
##  4.68074523 -0.08171004  0.38018558 
## 
##                Estimate      S.E.    CI95.lo   CI95.hi Ztest       Zval
## (Intercept)  4.68074523 0.2413345  4.2077383 5.1537521   *** 19.3952605
## x           -0.08171004 0.3087647 -0.6868778 0.5234577       -0.2646353
## y            0.38018558 0.3097641 -0.2269410 0.9873121        1.2273389
plot(ajus.INH.2,pause=FALSE,superimpose = FALSE)

plot(ajus.INH.2,pause=FALSE,how="persp", theta=-30,phi=40,d=4)

plot(ajus.INH.2,pause=FALSE,how="contour", theta=-30,phi=40,d=4)

pred2<-predict(ajus.INH.2)
plot(pred2)

# modelo cinhomogeneo on tendencia espacial NO lineal
require(splines)
## Loading required package: splines
#ajus.INH.3<-ppm(sim.INHa~bs(x,df=2)+bs(y,df=2),Poisson())

ajus.INH.3<-ppm(sim.INHa~poly(x,2)+poly(y,2),Poisson())

ajus.INH.3
## Nonstationary Poisson process
## 
## Log intensity:  ~poly(x, 2) + poly(y, 2)
## 
## Fitted trend coefficients:
## (Intercept) poly(x, 2)1 poly(x, 2)2 poly(y, 2)1 poly(y, 2)2 
##   4.6815325  -0.3972158  15.8736932   2.8181625  14.0666392 
## 
##               Estimate      S.E.   CI95.lo   CI95.hi Ztest       Zval
## (Intercept)  4.6815325 0.1019326  4.481748  4.881317   *** 45.9277186
## poly(x, 2)1 -0.3972158 2.5786078 -5.451194  4.656763       -0.1540427
## poly(x, 2)2 15.8736932 2.8529773 10.281960 21.465426   ***  5.5639046
## poly(y, 2)1  2.8181625 2.6344315 -2.345228  7.981553        1.0697422
## poly(y, 2)2 14.0666392 2.8831711  8.415728 19.717551   ***  4.8788777
plot(ajus.INH.3,pause=FALSE,superimpose = FALSE)

plot(ajus.INH.3,pause=FALSE,how="persp", theta=-30,phi=40,d=4)

plot(ajus.INH.3,pause=FALSE,how="contour", theta=-30,phi=40,d=4)

pred3<-predict(ajus.INH.3)
plot(pred3)

Como funciona KDE

library(manipulate)
x<-c(10,13,14)

GrafKDE<-function(h)
{
plot(density(x,bw=h,from=8,to=16),main=paste("Ventana =",h),ylab="Intensidad",ylim=c(0,0.8),lwd=3,xlab="")
points(x,c(0,0,0),pch=16,col="green")
n1<-function(x){dnorm(x,10,h)/3}
n2<-function(x){dnorm(x,13,h)/3}
n3<-function(x){dnorm(x,14,h)/3}
curve(n1, 7, 17, lwd=2, axes = FALSE, xlab = "", ylab = "",add=T,col="green",lty=3)
curve(n2, 7, 17, lwd=2, axes = FALSE, xlab = "", ylab = "",add=T,col="green",lty=3)
curve(n3, 7, 17, lwd=2, axes = FALSE, xlab = "", ylab = "",add=T,col="green",lty=3)
}
GrafKDE(0.5)
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "axes" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "axes" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "axes" is not a graphical
## parameter

#manipulate(GrafKDE(h=h),h=slider(0.2,2,step=0.1))

Aplicacion de KDE a Poisson Homogeneo

require(splancs)
## Loading required package: splancs
## 
## Spatial Point Pattern Analysis Code in S-Plus
##  
##  Version 2 - Spatial and Space-Time analysis
# eleccion de ventana por diggle
mserw <- bw.diggle(sim.CSR)
bw <- as.numeric(mserw)*3
#
den.CSR<-density(sim.CSR,bw)
plot(den.CSR)
points(sim.CSR)

Aplicacion de KDE a Poisson InHomogeneo

require(splancs)
# eleccion de ventana por diggle
mserw <- bw.diggle(sim.INHa)
bw <- as.numeric(mserw)*3
#
den.INH<-density(sim.INHa,bw)
plot(den.INH)
points(sim.INHa)

Elegiendo la ventana optima para Poisson Inhomogeneo

require(spatstat.utils)
## Loading required package: spatstat.utils
#datos<-sim.CSR
datos<-sim.INHa
poligono<-matrix(c(0,0,0,1,1,0,1,1,0,0),4,2,byrow = T)
## Warning in matrix(c(0, 0, 0, 1, 1, 0, 1, 1, 0, 0), 4, 2, byrow = T): data length
## [10] is not a sub-multiple or multiple of the number of rows [4]
Mse2d <- mse2d(as.points(datos), poligono, nsmse=50, range=1)
#png("/home/andresfaral/Dropbox/Estadistica Espacial/KDE INH CV")
plot(Mse2d$h,Mse2d$mse, type="l",xlab="Ventana (h)", ylab="Error Estimado")
arrows(0.15,9,0.15,-1.5,lty=1,col="green",lwd=3)

#dev.off()
bw<-min(Mse2d$h)
bw<-0.03
den.INH<-density(datos,bw)
png("/home/andresfaral/Dropbox/Estadistica Espacial/KDE INH chico")
plot(den.INH)
points(datos)
dev.off()
## png 
##   2

Elegiendo la ventana optima para Thomas

datos<-sim.THO
poligono<-matrix(c(0,0,0,1,1,0,1,1,0,0),4,2,byrow = T)
## Warning in matrix(c(0, 0, 0, 1, 1, 0, 1, 1, 0, 0), 4, 2, byrow = T): data length
## [10] is not a sub-multiple or multiple of the number of rows [4]
Mse2d <- mse2d(as.points(datos), poligono, nsmse=50, range=1)
bw<-Mse2d$h[which.min(Mse2d$mse)] 
png("/home/andresfaral/Dropbox/Estadistica Espacial/KDE THO CV")
plot(Mse2d$h,Mse2d$mse, type="l",xlab="Ventana (h)", ylab="Error Estimado")
arrows(bw,-1,bw,-4.5,lty=1,col="green",lwd=2)
dev.off()
## png 
##   2
den.THO<-density(datos,bw)
png("/home/andresfaral/Dropbox/Estadistica Espacial/KDE THO")
plot(den.THO)
points(datos)
dev.off()
## png 
##   2

Procesos Marcados

library(gstat)
## 
## Attaching package: 'gstat'
## The following object is masked from 'package:spatstat.core':
## 
##     idw
library(lattice)
## 
## Attaching package: 'lattice'
## The following object is masked from 'package:spatstat.core':
## 
##     panel.histogram
library(sp)
data(meuse)
class(meuse)
## [1] "data.frame"
coordinates(meuse) <- c("x", "y")
# EDA
spplot(meuse, "zinc", do.log = T, colorkey = TRUE)

bubble(meuse, "zinc", do.log = T, key.space = "bottom")

# fitting
xyplot(log(zinc) ~ sqrt(dist), as.data.frame(meuse))

zn.lm <- lm(log(zinc) ~ sqrt(dist), meuse)
meuse$fitted.s <- predict(zn.lm, meuse) - mean(predict(zn.lm,
meuse))
meuse$residuals <- residuals(zn.lm)
spplot(meuse, c("fitted.s", "residuals"))

Interpolacion por IDW

data(meuse.grid)
class(meuse.grid)
## [1] "data.frame"
coordinates(meuse.grid) <- c("x", "y")
meuse.grid <- as(meuse.grid, "SpatialPixelsDataFrame")
idw.out <- gstat::idw(zinc ~ 1, meuse, meuse.grid, idp = 0.5)
## [inverse distance weighted interpolation]
plot(idw.out)

Interpolacion con Modelado Lineal

#zn.lm <- lm(log(zinc) ~ sqrt(dist), meuse)
#zn.lm <- lm(zinc ~ sqrt(dist), meuse)
zn.lm <- lm(zinc ~ poly(x, y, degree = 2), meuse)
meuse.grid$pred <- predict(zn.lm, meuse.grid)
meuse.grid$se.fit <- predict(zn.lm, meuse.grid, se.fit = TRUE)$se.fit
plot(meuse.grid["pred"])

#image(meuse.grid["pred"])

El Variograma de zinc

hist(dist(coordinates(meuse)))

hscat(zinc ~ 1, meuse, seq(0,2000,length.out = 10))

#
sel <- plot(variogram(zinc ~ 1, meuse, cloud = FALSE))
sel

El Variograma de residuos

hscat(residuals ~ 1, meuse, seq(0,2000,length.out = 10))

#
sel <- plot(variogram(residuals ~ 1, meuse, cloud = FALSE))
sel

#
sel2<-variogram(log(zinc) ~ sqrt(dist), meuse)
plot(sel2)

#
zn.lm2 <- lm(zinc ~ dist+poly(x, y, degree = 2), meuse)
meuse$fitted2 <- predict(zn.lm2, meuse) - mean(predict(zn.lm2,
meuse))
meuse$residuals2 <- residuals(zn.lm2)
hscat(residuals2 ~ 1, meuse, seq(0,2000,length.out = 10))

#
sel <- plot(variogram(residuals2 ~ 1, meuse, cloud = FALSE))
sel

Procesos Marcados

SIN Autocorrelacion

require(gstat)
#
set.seed(1)
sim.CSR.M<-sim.CSR
X<-rnorm(sim.CSR.M$n)
X<-as.numeric(scale(X))
marks(sim.CSR.M)<-X
sim.CSR.M
## Marked planar point pattern: 91 points
## marks are numeric, of storage type  'double'
## window: rectangle = [0, 1] x [0, 1] units
sim.CSR.M.df<-as.data.frame(sim.CSR.M)
sim.CSR.M.spdf<-sim.CSR.M.df
coordinates(sim.CSR.M.spdf) <- c("x", "y") # como spatialpoints dataframe

plot(sim.CSR.M)

spplot(sim.CSR.M.spdf, colorkey = TRUE)

### Vaiograma
hscat(marks ~ 1, sim.CSR.M.spdf, seq(0,1,length.out = 11))

variog<-variogram(marks ~ 1, sim.CSR.M.spdf)
plot(variog,type="l")

CON Autocorrelacion

A proṕosito NO uso ni gstat ni fields

Genero valores al azar y luego hago muchas pasadas (pp) en las que elijo un punto al azar y modifico su marca sumandole una fracción (ff) de los valores de los (Kk) puntos cercanos

require(nabor) # Libreria para calcular vecinos mas cercanos 
## Loading required package: nabor
sim.CSR.M2<-sim.CSR.M
n<-sim.CSR.M2$n
pp<-500
ff<-0.5
kk<-3
#
coordenadas<-coords(sim.CSR.M2)
vecinos<-knn(coordenadas,coordenadas,k=kk+1)$nn.idx[,2:(kk+1)] # busca vecinos 
X.sac<-marks(sim.CSR.M2)
for (i in 1:pp)
{
  cual<-sample(1:n,1)
  X.sac[cual]<-X.sac[cual]+ff*mean(X.sac[vecinos[cual,]])
}
X.sac<-as.numeric(scale(X.sac))
marks(sim.CSR.M2)<-X.sac
sim.CSR.M2.df<-as.data.frame(sim.CSR.M2)
sim.CSR.M2.spdf<-sim.CSR.M2.df
coordinates(sim.CSR.M2.spdf) <- c("x", "y") # como spatialpoints dataframe

plot(sim.CSR.M2)

spplot(sim.CSR.M2.spdf, colorkey = TRUE)

### Vaiograma
hscat(marks ~ 1, sim.CSR.M2.spdf, seq(0,1,length.out = 11))

variog<-variogram(marks ~ 1, sim.CSR.M2.spdf)
plot(variog,type="l")

Indice de Moran

require(ape)
## Loading required package: ape
## 
## Attaching package: 'ape'
## The following object is masked from 'package:splancs':
## 
##     zoom
## The following objects are masked from 'package:spatstat.geom':
## 
##     edges, rotate
#
datos<-sim.CSR.M
datos.dists <- as.matrix(dist(cbind(sim.CSR.M$x, sim.CSR.M$y)))
datos.dists.inv <- 1/datos.dists
diag(datos.dists.inv) <- 0
Moran.I(marks(datos),datos.dists.inv)
## $observed
## [1] -0.04942335
## 
## $expected
## [1] -0.01111111
## 
## $sd
## [1] 0.01848804
## 
## $p.value
## [1] 0.03824008
#
datos<-sim.CSR.M2
datos.dists <- as.matrix(dist(cbind(sim.CSR.M$x, sim.CSR.M$y)))
datos.dists.inv <- 1/datos.dists
diag(datos.dists.inv) <- 0
Moran.I(marks(datos),datos.dists.inv)
## $observed
## [1] 0.1257112
## 
## $expected
## [1] -0.01111111
## 
## $sd
## [1] 0.01842769
## 
## $p.value
## [1] 1.130207e-13

Kriging

require(gstst)
## Loading required package: gstst
## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'gstst'
data(meuse.grid)
class(meuse)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
class(meuse.grid)
## [1] "data.frame"
coordinates(meuse.grid) <- c("x", "y")

kri.meuse<- krige(zinc~1, meuse, meuse.grid)
## [inverse distance weighted interpolation]
class(meuse.grid)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
plot(kri.meuse)

spplot(kri.meuse["var1.pred"])

#
M2.grid<-expand.grid(x=seq(0,1,length.out = 100),y=seq(0,1,length.out = 100))
coordinates(M2.grid) <- c("x", "y")
class(M2.grid)
## [1] "SpatialPoints"
## attr(,"package")
## [1] "sp"
m <- vgm(.59, "Sph", 874, .04)
kri.M2 <- krige(marks ~ 1, sim.CSR.M2.spdf,M2.grid,model=m)
## [using ordinary kriging]
plot(kri.M2)

plot(sim.CSR.M2)

spplot(kri.M2["var1.pred"], colorkey = TRUE)